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Abstract 

Stochastic reaction networks are dynamical models of biochemical reaction systems and form a 
particular class of continuous-time Markov chains on N”. Here we provide a fundamental character¬ 
isation that connects structural properties of a network to its dynamical features. Specifically, we 
define the notion of ‘stochastically complex balanced systems’ in terms of the network’s stationary 
distribution and provide a characterisation of stochastically complex balanced systems, parallel to 
that established in the 70-80ies for deterministic reaction networks. Additionally, we establish that 
a network is stochastically complex balanced if and only if an associated deterministic network is 
complex balanced (in the deterministic sense), thereby proving a strong link between the theory of 
stochastic and deterministic networks. Further, we prove a stochastic version of the ‘deficiency zero 
theorem’ and show that any (not only complex balanced) deficiency zero reaction network has a 
product-form Poisson-like stationary distribution on all irreducible components. Finally, we provide 
sufficient conditions for when a product-form Poisson-like distribution on a single (or all) compo¬ 
nent (s) implies the network is complex balanced, and explore the possibility to characterise complex 
balanced systems in terms of product-form Poisson-like stationary distributions. 


1 Introduction 

Improved experimental techniques have made it possible to measure molecular fluctuations at a small 
scale, creating a need for a stochastic description of molecular data [n na. Typically, biochemical 
reaction networks are modelled as deterministic systems of ordinary differential equations (ODEs), but 
these models assume the Individual species are in high concentrations and do not allow for stochastic 
fluctuation. An alternative is stochastic models based on continuous-time Markov chains [m nn m II 
Ena- As an example of a stochastic reaction system, consider 

A + B^2C, ( 1 . 1 ) 

K2 

where Ki,K 2 are positive reaction constants. The network consists of three chemical species A, B and 
C and two reactions. Each occurrence of a reaction modifies the species counts, for example, when the 
reaction A-\- B ^ 2C takes places, the amount of A and B molecules are each decreased by one, while 
two molecules of C are created. The species counts are modelled as a continuous-time Markov chain, 
where the transitions are single occurrences of reactions with transition rates 

Xl{x) = KiXAXB, X2{x) = K2Xcixc - 

and X = {xa,xb,xc) are the species counts [1]. When modelled deterministically, the concentrations 
(rather than the counts) of the species change according to an ODE system. 

In a classical paper m, Kurtz explored the relationship between deterministic and stochastic reaction 
systems, using a scaling argument - large volume limit - to link the dynamical behaviour of the two 
types of systems to each other. Other, mainly recent work, also points to close connections between the 
two types of systems miiiiiiisiiis]. In this paper we explore this relationship further. 
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by the Danish Research Councils. 
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A fundamental link between structural network properties and dynamical features of deterministic 
reaction networks has been known since the 1970s and 1980s with the work of Horn, Jackson and 
Feinberg [laiio]. Specifically, their theory concerns the existence and uniqueness of equilibria in complex 
balanced systems, with the ‘deficiency zero theorem’ playing a central role in this context. Complex 
balanced systems were called cyclic balanced systems by Boltzmann. They have attractive analytical 
and physical properties; for example a (pseudo-)entropy might be defined which increases along all 
trajectories (Boltzmann’s H-theorem) [71 lid). 

A parallel theory for the stochastic regime is not available, and the very concept of “complex balanced” 
does not currently have a stochastic counterpart. In this paper we develop a theory to fill this gap. We 
define stochastically complex balanced systems through properties of the stationary distribution, and 
we prove results for stochastic reaction networks that are in direct correspondence with the results for 
deterministic models. In particular, we prove a parallel statement of the deficiency zero theorem and 
show that all deficiency zero reaction networks have product-form Poisson-like stationary distributions, 
irrespectively whether they are complexed balanced or not. In fact, in the non-complexed balanced case, 
the network is complex balanced on the boundary of the state space. 

A second target of our study concerns product-form stationary distributions. Such distributions 
are computationally and analytically tractable and appear in many areas of applied probability, such 
as, queueing theory [13 im, Petri Net theory m, and stochastic reaction network theory [281 1201 
[2]. Specifically, a complex balanced mass-action network has a product-form Poisson-like stationary 
distribution on every irreducible component unii]. As an example, the stationary distribution of the 
complex balanced reaction system is 

^XA^XB ^XC 

TTr{x) = Mr - - ■ - ■ forager, 

xa'-xb'-xc'- 

where P = {a: S : a;^ + a;_B -I- 2a;c = 0} is an irreducible component of the state space and Mr is a 
normalising constant. 

We expand the above result on mass-action systems and give general conditions under which the 
converse statement is true. In particular, we are interested in providing a structural characterisation of 
the networks with product-form Poisson-like stationary distributions. However, this class of networks is 
strictly larger than that of complex balanced networks, and a full characterisation seems hard to achieve. 
We illustrate this with examples. 

2 Background 

We first introduce the necessary notation and background material; see mdQiE] for general references. 
We assume standard knowledge about continuous-time Markov chains. 

2.1 Notation 

We let M, Kq and K+ be the real, the non-negative real and the positive real numbers, respectively. Also 
let N be the natural numbers including 0. 

For any real number a G K, |a| denotes the absolute value of a. Moreover, for any vector v G 
we let Vi be the ith component of v, ||u|| the Euclidean norm, and ||u||oo the infinity norm, that is, 
||u||oo = maxi |ui|. For two vectors v,w G K^, we write v < w (resp. u > w) and v < w (resp. u > w), if 
the inequality holds component-wise. Further, we define l{„<uj} to be one if u < w, and zero otherwise, 
and similarly for the other inequalities. If u > 0 then v is said to be positive. Finally, supp v denotes the 
index set of the non-zero components. For example, if u = (0,1,1) then suppu = {2, 3}. 

If X G ffig and u G we define 
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with the conventions that 0! = I and 0° = I. 
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2.2 Reaction networks 


A reaction network is a triple G = {X,C, TZ), where X = {Si,S 2 , ■ ■ ■, Sn} is a set of n species, C is a set 
of m complexes, and 7?. C C x C is a set of fc reactions, such that (y, y) ^ TZ for all y GC. The complexes 
are linear combinations of species on N, identified as vectors in K". A reaction {y,y') G TZ is denoted by 
y ^ y'. We require that every species is part of at least one complex, and that every complex is part of 
at least one reaction. In this way, there are no “superfluous” species or complexes and G is completely 
determined by the set of reactions 7Z, which we allow to be empty. In there are n = 3 species 

{A, B,C), m = 2 complexes {A + B, 2C), and k = 2 reactions. 

Given a reaction network G^ the reaction graph of G is the directed graph with node set C and edge set 
TZ. We let i be the number of linkage classes (connected components) of the reaction graph. A reaction 
y ^ y' G TZ is terminal if any directed path that starts with y —> y' is contained in a closed directed 
path. We let TZ* be the set of terminal reactions. 


A reaction network G is weakly reversible, if every reaction is terminal. The network in (1.1) is weakly 
reversible, since both reactions are terminal. 

The stoichiometric subspace of G is the linear subspace of K" given by 

S = span(y' - y|y ^ y' G TZ). 

For V G K", the sets (?; + S') n Kq are called the stoichiometric compatibility classes of G (Fig. lA). For 
the network in (1.1), S = span((—1,1, 0), (0,1, —1)) C which is 2-dimensional. 


2.3 Dynamical systems 

We will consider a reaction network G either as a deterministic dynamical system on the continuous 
space Kq , or as a stochastic dynamical system on the discrete space N". 

In the deterministic case, the evolution of the species concentrations z = z{t) G Rq at time t is 
modelled as the solution to the ODE 

^ (2/'-2/)Ay^y'(z), (2.1) 

for some functions Xy^yi : Kq — )■ Kq and an initial condition z(0) G Kq . We require that the functions 
Xy^y' are continuously differentiable, and that \y^yi{z) > 0 if and only if suppy C suppz. Such 
functions are called rate functions, they constitute a deterministic kinetics K for G, and the pair {G,K) 
is called a deterministic reaction system. If Xy^yi{z) = Ky^y'Z^ for all reactions, then the constants 
Ky^y' are referred to as rate constants and the modelling regime is referred to as deterministic mass- 
action kinetics. In this case, the pair {G, k) is called a deterministic mass-action system, where k G 
is the vector of rate constants. 

In the stochastic setting, the evolution of the species counts X(t) G N" at time t is modelled as a 
continuous-time Markov chain with state space N". At any state x G N", the states that can be reached 
in one step are x -\- y' — y tor y ^ y' G TZ, with transition rates Xy^y' (x) . The functions Xy^y /: N" — >■ Kq 
are called rate functions, and we require that Xy^yi{x) > 0 if and only if a; > y. A choice of these 
functions constitute a stochastic kinetics K for G and the pair (G, K) is called a stochastic reaction 
system. If the reaction y ^ y' occurs at time t, then the new state is 

X{t) = X{t-) -\-y' -y, 


where X{t—) denotes the previous state. If for any reaction y ^ y' GTZ 

T.\ 

then the constants Ky^yi are known as rate constants, as in the deterministic case, and the modelling 
regime is referred to as stochastic mass-action kinetics. The pair (G, k) is, in this case, called a stochastic 
mass-action system. 

The evolution of the stochastic as well as the deterministic reaction system is confined to the stoi¬ 
chiometric compatibility classes. 


A 


y^y' 


(x) = 


v^y' 


(x- 


z(t) G (z(0) -f S') n and A(t) G (A(0) -f S) n 
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In fact, X{t) G (-^(0) + S') n N", as X{t) takes values in N”. 

Definition 1. Let Q = {X,C,TZ) be a reaction network. 

a) A reaction network Q' = {X',C',TZ') is a subnetwork of Q if TZ' Q TZ. In this case, it follows that 
X' CX andC CC. 

b) A system {Q', K'), deterministic or stochastic, is a subsystem of a system {Q, K) ifQ' is a subnetwork 
of Q and the rate functions agree on the reactions in TZ'. 

c) The subnetwork Q* given by the set of terminal reactions TZ* is the terminal network ofQ. We denote 
Q* = (X*,C*,TZ*). Furthermore, the subsystem {Q*,K*) of {Q,K) is called the terminal system of 
{G,K). 

Definition 2. The connected components of the reaction graph of the terminal network of Q are called 
terminal strongly connected component of Q. For any complex y in C*, we denote by (Qy,Ky) the 
subsystem of Q whose reaction graph is the terminal strongly connected component containing y as node. 

As an example, consider the mass-action system 

2A^2B<^ AZ^O^C. 

K2 K-e 

Here, there are two terminal strongly connected components, which are 2A ^ 2B and 0 ^ C. In 
particular, {G 2 At K 2 a) is equal to {G 2 B,K 2 b) and is given by 

2A ^ 2B. 

K2 

Finally, if (G,k) is a mass-action system, any subsystems {G',K') is a mass-action systems as well and 
can be denoted by {G', k'). 


3 Deterministic reaction systems 

In this section we will recapitulate the known characterisation of existence and uniqueness of positive 
equilibria in complex balanced systems and the connection between complex balanced systems and 
deficiency zero reaction networks. As we will show in the subsequent section, this characterisation can 
be fully translated into a similar characterisation for stochastic reaction networks. 


3.1 Complex balanced systems 

We start with a definition. 


Definition 3. A deterministic reaction system (G,K) is said to be complex balanced if there exists a 
positive complex balanced equilibrium, that is, a positive equilibrium point c G K" for the system (2.11, 
such that 

^ Xy^y-ic) = ^ Xy'^yic) foc all y € C. (3.1) 


y'ec 


y'ec 


The name ‘complex balanced’ refers to the fact that the flow, at equilibrium, entering into the complex 


y equals the flow exiting from the complex. As an example, the mass-action system in (I.I) is complex 


balanced for any choice of (^ 1 ,^ 2 ) and c = (k 2 ,«^i,«;i) is a complex balanced equilibrium. The class of 
complex balanced systems is an extension of the class of detailed balanced mass-action systems [muD]. 
For mass-action systems, (|3.1[) becomes 


E' 

y'ec 


'‘V^y 


^y _ 


E 

v'ec 


Kyi^yC^ for all y G C, 


(3.2) 


with the convention that ky^yi = 0 ii y ^ y' ^ TZ. 

In the case of mass-action kinetics, we extend Definition to the stochastic case, by saying that a 
stochastic mass-action system {G, k) is complex balanced if the deterministic mass-action system (G, k) is 
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complex balanced. We might therefore refer to complex balanced mass-action systems without specifying 
whether they are stochastically or deterministically modelled. 

The next theorem is a slight generalization of a classical result m, which provides the backbone for 
the further characterisation. The generalization includes a property of non-negative equilibria 


Theorem 3.1. If a deterministic reaction system {Q, K) is complex balanced, then Q is weakly reversible. 
Moreover, if K is mass-action kinetics, then all equilibria are complex balanced, that is, they fulfil (3.2|. 
Moreover, there exists exactly one positive equilibrium in each stoichiometric compatibility class, which 
is locally asymptotically stable. 

As we are not aware of a proof of this more general formulation, we provide one in Appendix [B| 


3.2 Deficiency zero statements 

The deficiency plays an important role in the study of complex balanced systems. The dehciency of Q 
is defined as 

S = m — £ — s, 

where m is the cardinality of C, £ is the number of linkage classes of the reaction graph of Q and s is the 
dimension of the stoichiometric subspace S |13j . The definition hides the geometrical interpretation of 
the deficiency, which we now will explore. 

Let {cylygc be a basis of K™. Further, define 

dy — ‘,y' — Cyf Cy and ^y—^yf — y y 


for y ^ y' G TZ. Let D = span{dy^y'\y —>■?/'€ TZ). Then dimH = m — £ [13]. 

The space D is linearly isomorphic to the stoichiometric subspace S if and only if 5 = 0. Specifically, 
consider the homomorphism 

(p: M™ -)► K” 

Cy ^ y. 

For y ^ y' G TZ, we have (p{dy^y') = q:>\D- D ^ S is thus a surjective homomorphism. 

Therefore, 

dim Ker <p\d = dim D — s = m — £ — s = 6, (3.4) 

which implies that is an isomorphism if and only if 5 = 0. It further follows that the deficiency is a 
non-negative number. 

We state here a useful Lemma on the deficiency of subnetworks. 


Lemma 3.2. Let Q be a reaction network with deficiency 6. Then, the deficiency of any subnetwork of 
Q is smaller than or equal to 6. 


Proof. Let TZ' <GTZ and let Q' be the corresponding subnetwork with deficiency 6'. Further, let D' and 
S' be the equivalent of D and S for Q', respectively. By (3.4) and since D' is a subspace of D, we have 
S' = dimKer(/j|£)/ < dim Ker = S, which concludes the proof. □ 


We next state two classical results which elucidate the connection between complex balanced systems 
and deficiency zero systems. A proof of the first and of the second result can be found in m and in m, 
respectively. The results draw a connection between graphical and dynamical properties of a network. 
Theorem 3.4 is given here in a wider formulation than in [lOj (see Appendixfor a proof). 


Theorem 3.3. The mass-action system {Q, k) is complex balanced for any choice of k if and only if Q 
is weakly reversible and its deficiency is zero. 


Theorem 3.4. Consider a deterministic reaction system {Q,K), and assume that the deficiency of Q is 
zero. If X G Kq is an equilibrium point and y ^ y' GiZ, then supp y C suppx only ify^y' is terminal. 
Moreover, if K is mass-action kinetics with rate constants k and suppy C suppx, then the projection of 
X onto the species space of Qy is a complex balanced equilibrium of {Gy, Ky). 


It follows from Theorem 3.4 that an equilibrium point satisfies (3.2) for the terminal system, though 
it is not necessarily a positive equilibrium of (t/*,K*).The deficiency zero theorem, in the following 
formulation, is a consequence of the three previous theorems: 
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Figure 1: The figure shows some features of the reaction network 2A —>■ 2B, and A + 35 —>■ 3A + B. 
(A) The stoichiometric compatibility classes are of the form {{za,zb)' za + Zb = const.}. (B) The 
two irreducible components on {{xa^xb)- Xa + Xb =6} are shown (black circles and square), together 
with the possible transitions between the states. All states within a component are accessible from each 
other. The “square” component has no active reactions, both reactions are active on the “black circles” 
component. The grey states are transient states which are not in any irreducible component. 


Theorem 3.5 (Deficiency zero theorem). Consider a deterministic reaction system {Q,K) for which the 
deficiency is zero. Then the following statements hold: 

i) if G is not weakly reversible, then there exists no positive equilibria; 

ii) if G is weakly reversible and K is mass-action kinetics, then there exists within each stoichiometric 
compatibility class a unique positive equilibrium, which is asymptotically stable. 

The original formulation is richer than the one presented here uni- 

4 Stochastic reaction systems 

4.1 Classification of states and sets 

To characterise the stochastic dynamics we introduce the following terminology. 

Definition 4. Let G = {X,C,TZ) be a reaction network. 

a) A reaction y ^ y' G TZ is active on x G N” if x>y. 

b) A state u G N” is accessible from a state x G N" if there is a sequence of q > 0 reactions (pj -G 
y'j)j=i,...,q such that 


(i) u = x + -%■); 


(ii) yh -t y( is active on x-G X]j=i “ %) all 1 < h < q. 

Definition 5. Let G be a reaction network. A non-empty set T C N” is an irreducible component of G 
if for all X GT and all u G N", u is accessible from x if and only if u GT. 


Definition 6. A reaction network G is essential if the state space is a union of irreducible components. 
A reaction network G is almost essential if the state space is a union of irreducible components except 
for a finite number of states. 
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An essential network is also almost essential. A weakly reversible reaction network is essential [52]. 
Conditions for being essential can be found in [221 Ej. Any irreducible component is contained in 
some stoichiometric compatibility class, and a stoichiometric compatibility class may contain several 
irreducible components (Fig. IB). 

4.2 Stationary distribution 

The stationary distribution 7 rr on an irreducible component F is unique, if it exists. It is characterised 
by the master equation [4]: 

Trr{x + y-y')Xy^y>{x + y-y')=TTr{x) Y (4-1) 

y-i-y'GT?, y-j-i/'GT?, 

for all a; G F. Let X{t) denote the stochastic process associated with the system. If X{to) follows the 
law of TTp at time t^, then the distribution of X(t) is TTp for all future times t > t^. In this sense, the 
stationary distribution describes a state of equilibrium of the system. Moreover, if 7 rr exists, then 

lim P{X{t) G A) = 7 rr(A) for any A C F, (4-2) 

t—^OO 

provided that A(0) G F with probability one. As discussed in Section a connection between mass- 
action complex balanced systems and their stationary distribution has been made in [ 2 ]: 

Theorem 4.1. Let {G,k) be a complex balanced mass-action system. Then, there exists a unique sta¬ 
tionary distribution on every irreducible component F, and it is of the form 

ttt{x) = /o’’ 2 ;GF, (4.3) 

where c is a positive complex balanced equilibrium of {Q, n) and Mf is a normalising constant. 

4.3 Parallel theorems for stochastic mass-action systems 

In this section we derive stochastic statements corresponding to Theorem |3.1|3.5| Some of the proofs 
are deferred to Appendix |B| We begin with a definition. 

Definition 7. For an irreducible component F, the set TZr of active reactions on F consists of the 
reactions y —>■ y' G TZ that are active on some a; G F. The subnetwork Qr = {Xy',Ct,'R-t) is called the 
F-network of G and the subsystem (Gr^Kr) of{G,K) is called the F-system of{G,K). 

The reactions that are active on F determine the dynamics of the stochastic system on F. To study 
the stationary distributions, it is therefore convenient to analyse the F-systems. Note that TZy is empty 
if and only if F consists of a single state. 

As an example, consider the deficiency zero network, 

C ^ D, 2A^2B, A ^ 0. 

All molecules of A and B are irreversibly consumed through A —>■ 0 and 2B —)■ 2A, thus the only active 
reactions on an irreducible component F 7 ^ {0} are C ^ D. The T-network is therefore C ^ D, which 
differs from the terminal system C ^ D, 2A ^ 2B. The next proposition states that for a deficiency 
zero reaction network T^r ^ TZ* for any irreducible component F. Note that Proposition |4.2| does not 
hold in general, for example, 

A^ B, 2B ^2A 

has TZy = TZ for any F 7 ^ {0}, {(0,1)}, while TZ* = 0. 

Proposition 4.2. Let G be a reaction network and F an irreducible component such that Gr has deficiency 
zero. Then, Gr is a subnetwork of G* ■ In particular, this is true if the deficiency of G is zero. 

See Appendix [^ for a proof. Proposition |4.2| can be useful because TZr might be difficult to find, 
especially if there are many complexes. On the other hand, terminal reactions are easily identified by 
means of the reaction graph. The next definitions are inspired by Definition 
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Definition 8. Let {Q, K) be a stochastic reaction system. A stationary distribution nr on an irreducible 
component T is said to be complex balanced if 

nr{x - y'+y)Xy^y>{x - y'+ y) = '^ nr{x)Xy'^y{x) Vy'eCr,a;er. (4.4) 

veCr yGCr 


For a mass-action system, (4.4) becomes 


^ nr{x-y' + y) 
yeCr 


(x-y' + y)l 

(x-yy. 


xl 


veCr y >■ 


for any y' S Cr and a: G F, with the convention that ky^yi = 0 ii y ^ y' ^ TZr- In developing 
the theory for complex balanced equilibria in the deterministic setting, an important role is played by 
requiring positivity of the complex balanced equilibrium. Our aim is to introduce a similar concept for 
the stochastic systems. In the deterministic setting, if a state z S M" is positive then every rate function 
calculated on z is positive. We find inspiration from this to give the next definition: 

Definition 9. An irreducible component F is positive ifQr = G- 

Equivalently, an irreducible component F is positive if all reactions are active on F. The next definition 
follows naturally by analogy with the deterministic setting. 


Definition 10. A stochastic reaction system (G, K) is said to be stochastically complex balanced if there 
exists a complex balanced stationary distribution on a positive irreducible component. 



Theorem 4.3. Let {G,K) be a stochastic reaction system, and let F be an irreducible component. If 
there exists a complex balanced stationary distribution nr on F then Gr is weakly reversible. Moreover, if 
K is mass-action kinetics with rate constants k, there exists a complex balanced stationary distribution 
nr on F if and only if the T-system of (f/, k) is complex balanced. If this is the case, then nr has the 
form 

_ Xi 

nr{x) = Mf ^ /orxSF, (4.5) 

where c is a positive complex balanced equilibrium of {Gr, Kr) and Mf is a normalising constant. 

The proof is in Appendix]^ ft is shown in [5] that the stationary distribution 7rr(x) is independent 
of the choice of complex balanced equilibrium c of the F-system, provided that it is positive. We are now 
ready to derive stochastic versions of Theorem |3.1||3.5| In addition, we will show that a stochastically 
complexed balanced mass-action system is complex balanced and vice versa. Hence, we will show that 
the deterministic and stochastic systems are intimately connected. The next corollary is an analogue of 
Theorem 13.11 


Corollary 4.4. If a stochastic reaction system {G, K) is stochastically complex balanced then G is weakly 
reversible. Moreover, a mass-action system {G, k) is stochastically complex balanced if and only if it is 
complex balanced. If this is case, then on every irreducible component F there exists a unique stationary 


distribution nr- Such 7rr is a complex balanced stationary distribution and it has the form (4.3), where 
c is a positive complex balanced equilibrium of {G,k). 


Proof. If F is positive, then (tJrjATr) = {G,K)- Therefore, by Theorem 4.3 if {G,K) is stochastically 
complex balanced then G is weakly reversible. Moreover, if K is mass-action kinetics with rate constants 
K, it follows from Theorem |4.3| that there exists a complex balanced stationary distribution on F if and 


only if {G, n) is complex balanced. In this case, by Theorem 4.1 a stationary distribution exists on every 
irreducible component and it is of the form (4.3). By Theorem 4.3 it is a complex balanced stationary 
distribution. □ 
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Corollary 4.4 might be considered a stochastic version of Theorem 3.1 especially if (4.21 is taken to 
be equivalent to “asymptotic stability” for a deterministic equilibrium. Part of the corollary is known 
[2] (see also Theorem 4.1), and the whole corollary might therefore be considered as an extension of the 
result in [5] on mass-action systems. In this sense, Theorem 4.3 provides an even more general version, 
which deals with complex balanced subsystems of {Q, k). 

We now state the parallel versions of Theorem |3.3||3.5| for the stochastic setting. 


Corollary 4.5. The mass-action system {G,k) is stochastically complex balanced for any choice of k, if 
and only if G is weakly reversible and its deficiency is zero. 


Proof. The result is an immediate consequence of Corollary |4.4| and Theorem |3.3| □ 

Theorem 4.6. Consider a stochastic reaction system {G,K), and assume the deficiency of G is zero. 
Let X be a state in an irreducible component T and let y ^ y' in TZ. Then, y < x only if y ^ y' is 
terminal. Moreover, if K is mass-action kinetics, then on T the stationary distribution has the form 


TTrix) = Mf 

i: SiGX’^ 


Xi\ 


for a: S r, 


(4.6) 


where c is a positive complex balanced equilibrium for the terminal system, and Mf is a normalising 
constant. 


The proof is in Appendix [B} 

Theorem 4.7. Consider a stochastic reaction system {G,K), and assume that the deficiency of G is 
zero. Then the following statements hold: 


i) if G is not weakly reversible, then there exist no positive irreducible components; 


ii) if G is weakly reversible, then G is essential, and if K is mass-action kinetics then there exists a 
unique stationary distribution on every irreducible component. 


The proof of the theorem is in Appendix]^ In case (i). Theorem 4.6 provides the form of the stationary 
distribution. Hence we have characterised the stationary distribution for any deficiency zero reaction 
system, irrespectively whether it is complex balanced or not. 


Example 1. Consider the two stochastic mass-action systems 



K2 


lOA ^ lOH and B, 


lOA ^ 0. 


The behaviours of the two corresponding deterministic systems differ substantially, while the behaviours 
of the stochastic systems are equivalent on the irreducible components Tg = {x S : xi -I- ^2 = 0} with 
0 < 0 < 10 an integer. Indeed, in both cases the Tg-system is 



K2 


which is 
on Tg is 


complex balanced (Theorem |3.3[). It follows from Theorem 4.3 that the stationary distribution 


Xi 

n,{x^,X2) = Me^^ 

Xi\ X2\ 


for (a;i,a; 2 ) S Pe, 


for a suitable normalizing constant Mg. The stationary distributions are complex balanced, but since 
Pg is not positive in either of the two networks, we cannot conclude that the systems are stochastically 
complex balanced. Indeed, they are not for some choice of rate constants (Corollary |4.5[ ). 

Incidentally, note that the second network is not almost essential. 
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5 Product-form Poisson-like stationary distributions 

The above results draw parallels between stochastic and deterministic reaction networks. If a mass- 
action system is (stochastically) complex balanced, then the stationary distribution on every irreducible 
component is a product-form Poisson-like distribution. Does the reverse statement hold true too? If the 
stationary distribution is a product-form Poisson-like distribution on some, or all irreducible components, 
does it follow that the system is complex balanced? In the spirit of the first part of the paper we would 
like to achieve a full characterisation of stochastic systems with product-form Poisson-like stationary 
distributions. However, even though the hypothesis of Theorem 5.1 below is rather general, a full 
characterisation seems hard to achieve. 

Theorem 5.1. Let Q be an almost essential reaction network, k G a vector of rate constants and 
c G K" a vector with positive entries. The probability distribution TTp: T —>■ (0,1], defined by (4.3) is a 
stationary distribution for the stochastic mass-action system (Q, k) for all irreducible components T C N" 
of Q if and only if c is a complex balanced equilibrium for (Q,k). 


Proof. By Theorem 4.1 if c > 0 is a complex balanced equilibrium for {Q,k), then the stationary 
distribution on all irreducible components P C N” is of the form (4.3). 

Oppositely, assume that (4.3) is the stationary distribution on P for the stochastic mass-action system 
{Q, k), for all irreducible components P. Since Q is almost essential, there exists a constant K such that 
any states x with ||a:|| > K belongs to an irreducible component P. For any x € N", such that 


mm Xi > max 


+ II2/II00) + K, 


(5.1) 


we have that x > y and x — y' y > y ioi all y ^ y' G TZ. Then, since (4.3) is a stationary distribution 
and since x and x -\- y — y' are in the same irreducible component for all y ^ y' G TZ, we have from (4.1) 

{x-\-y- y')\ 


T^vix + y-y')!- 

v^v'en 


{x - y')! 


= 7rr(x) 


y^y’eTl 


(x- y)!’ 


(5.2) 


for all X S P satisfying ( |5.1[ ). Further, using (4.3), equation ( |5.2[ ) becomes 

r* I 


E 

y-fy'^TZ 


Xl 


ix-y')n^^ 


'C 


y-y 


'= E 


x\ 




y-fy'^TZ 


(x-y)!’ 


which, by rearranging terms, leads to 

x! 


^ (x — y')\ ^ 


y'ec 


y^y'en 


^y^y'C - 2 ^ ( 

y'ec^^ y )■ y'^yen 


v'^y 


(5.3) 


The equality holds for all x S N" satisfying (5.1), therefore the polynomials on the two sides of (5.3) are 
equal. 


For any y' G C, let Py'{x) be the polynomial 

Py' (a:) = 


x\ 

(x - y')!' 


The monomial with maximal degree in py/ is x^ , and these differ for all complexes y' G C. This implies 
that Py!, y' G C, are linearly independent on K, and thus, the polynomials on the two sides of (5.3) are 
equal if and only if 


E' 


^y^y 


r^y-y 


' = E' 


“y ^y 


for all y' gC. 


yec yee 

Hence, c is a complex balanced equilibrium for {Q, k) and the proof is completed. 


□ 
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5.1 Relaxation of Assumptions in Theorem 5.1 


To infer the existence of positive complex balanced equilibria in Theorem |5.1[ the assumptions of the 


theorem could be weakened. Specifically, it is only required that (5.31 holds for a set of states whose 


geometry and cardinality allow us to conclude that the polynomials on the two sides of (5.31 are the same. 


For (5.3) to hold, we need x to be in a irreducible component and we require x > y and x — y' + y > y 
for all reactions y ^ y' G TZ, as well as the stationary distribution evaluated in x and x — y' + y to he oi 


the form (4.3). If a state x satisfies this, we call it a good state. 


A more general condition than being almost essential could be chosen case by case and depends on 


the monomials appearing in (5.3). For example, if the set of complexes coincides with the set of species. 


then the polynomials in (5.3) are linear and the existence of n + 1 good states in general position implies 


the existence of a positive complex balanced equilibrium. In general, let d be the total degree of the 
polynomials in (5.3). Then it is sufficient to have n lines in general position with more than d + 1 good 


states on each of them. Therefore, to conclude that a system is complex balanced it is sufficient to check 
the behaviour of a finite number of states, lying on a hnite number of irreducible components. However, 
it follows from Examples and that the existence of arbitrarily many good states on a few irreducible 
components does not imply the existence of a positive complex balanced equilibrium in general. Finally, 
in order to postulate that the mass-action system is complex balanced, it is necessary that the vector c 
appearing in Theorem |5.1| is the same for every irreducible component, as shown in Example]^ 

The following examples are also meant to give an idea of why it is hard to obtain a full characteriza¬ 
tion of stochastic mass-action systems with a product-form Poisson-like stationary distribution on some 
irreducible component. 


Example 2. Let p G K+ and let 6* > 2 be an integer. Consider the stochastic mass-action system 


A^B 


2B 


2A, 


(5.4) 


where ki = p{6 — 1) and K 2 = p are the rate constants. The reaction network is almost essential. 
It is shown in Appendix that the stationary distribution on the irreducible component Fg = {a: € 
: a;i -I- a ;2 = 6} has the form (4.3) with c = (1,1), namely 


T^e{xi,X 2 ) = Mg 


1 


a;i!a;2! 


for (a;i,a:: 2 ) G Fg, 


(5.5) 


where Mg is a normalising constant. However, the mass-action system is not complex balanced as the 
reaction network is not weakly reversible (Theorem |3.1[ ). In particular, by Theorem 5.1 not all irreducible 
components can have a stationary distribution of the form (4.3) with c = (1,1). Trivially, the absorbing 
states (0, 0) and (0,1) have it. 

Additionally, we should point out that there is not an equivalent system on F (that is, a stochastic 
mass-action system with the same transition rate matrix on the states of F as (|5.4[)) which is complex 


balanced. Consider the case 0 = 1. Since the transition from (0, 2) to (2,0) is possible according to (5.4), 
any equivalent mass-action system must contain the reaction 2B -G 2A, with rate constant p. It can 
be further shown that any equivalent weakly reversible mass-action system must contain the connected 
component 


A + B 



2A . 


This prevents the system from being complex balanced, since there is not a c G 7Z^ fulfilling (3.2) for 
the three complexes 2B, 2A and A + B. 


Example 3. Let pi,p 2 ,ps G M+ and let 0 > 2 be an integer. Consider the modihcation of Example 
given by 


A 


Pl(0-1) + P2 


B 


2B 


P1+P3 


2A, 


P2 


P3 


which is weakly reversible. If we let p 2 = 0 and ps = 0, then the system reduces to that of Example 
by removing the two reversible reactions. It can be shown that for any parameter choice, (|5.5[) is still 
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a stationary distribution on the irreducible component Fg = {x G : xi + X 2 = &}■ However, for 
some choice of parameters the mass-action system is not complex balanced. This can be seen either by 


direct computation on the system of complex balance equations (3.21 or by noting that the deficiency 


of the network is 1, so there must be a choice of parameters which prevents positive complex balanced 
equilibria by Theorem |3.3| It can be further shown that irreducible components different from T do not 
possess a product-form Poisson-like stationary distribution. 

Example 4. Consider the stochastic mass-action system with p G M+ and 9i, 02 two positive integers, 


A 

3H 


pQ\Q‘2 


^ B 


2B 2yl 


-G A+ 2B 2A + B 


-G 3B. 


The reaction network is almost essential. For any 0 G N, consider the irreducible component Fg = {cc G 
: xi -I- X 2 = 0 -I- 1}. Then TTg^ and irg^, defined as in (5.5), are the (unique) stationary distributions 
on the irreducible components Fgj^ and Fg^, respectively. For the relevant calculations see Appendix 
0 However, the mass-action system is not complex balanced, since the reaction network is not weakly 
reversible (Theorem |3.1[). 


Example 5. Theorem |5.1| can be also used to compute the stationary distribution of a stochastic mass- 
action system which behaves as a complex balanced system on the irreducible components. Consider the 
weakly reversible (and therefore essential) stochastic mass-action system 


A^2A 


A + B ^2A + B. 

K4 


On every irreducible component Fg = {x G : X2 = 0}, 0 G N, the associated continuous time 
Markov chain, which describes the evolution of the counts of A, has the same distribution as the process 
associated with 


K,2+K,i6 


because the transition rates coincide. The latter system is complex balanced for any choice of rate 
constants. The stationary distribution has the form (Theorem |5.1[) 


-Ke{x) = Mg — 


K2 + K 40 Y 
Kl + K3O ) 


for some positive constant Mg. The latter gives the stationary distribution of the original system as well. 
However, the rate of the Poisson distribution does depend on 9, in which case the original system cannot 
be complex balanced (Corollary |4.4[). For the same reason the example does not contradict Theorem 5.1 


6 Applications 

There are not many means to explicitly calculate the stationary distribution of a stochastic mass-action 
system. As an example, Theorem |4.3| can be used to determine the stationary distributions of mass-action 
systems like 

D, 2A^2B, A ^ 0. 

K 2 K4 

Indeed, for any irreducible component F different from {0}, the F-system is given by 



K2 


D, 


which is weakly reversible and has deficiency zero, therefore it is complex balanced. Hence, the stationary 
distribution on F has the form 


7rr(a;) = Mp 


0 : 3 ! X 4 ! 


for X G F, 


12 












where X 3 and X 4 denote the entries relative to C and D, respectively. Alternatively, since the terminal 
system is given by 




2A ^ 2B, 


K2 «4 

Theorem |4.6| can be used to compute the stationary distribution. On every irreducible component T, it 
is given by 


TTr{x) = Mr- 


'-1 


for a; € r, 


Xi\ X 2 \ Xz\ 2 : 4 ! 

which is equivalent to the previous formula since xi and X 2 are constantly 0 on all irreducible components. 
If the system does not fulfil the conditions of Theorem |4.3| and neither can be cast as a birth-death 


process, Theorem 5.1 might be useful. The following mass-action system is considered in [T]: 

A^O 0^2A. 


By Theorem 5.1 the stationary distribution cannot be Poisson. Indeed, it is given by the distribution 


of F = Yi -I- 2 Y 2 , where Fi and Y 2 are two independent Poisson random variables with rates ^ and 
respectively. Hence, 


7r(x) = e-^ y 

z!j! / \ 2 ki 




In the following system is also considered: 


O^A 

1^2 


2A ^ 3A. 


It has the stationary distribution 


r(x) =mA 


2=1 


ei[{i-i){i-2) + e2] 

i{i — l)(i — 2) -I- Bsi 


for X € N, 


where 9i = K^jK^, 02 = Kifn^, 9^ = K 2 /K 4 and M = 7 r( 0 ) is a normalising constant. It is interesting 
that 7 r(x) is a Poisson distribution if and only if 02 = 03 . In fact, and in accordance with our results, the 
mass-action system is complex balanced if and only if 02 = 03 - 

7 Discussion 


Corollary |4.5| provides a characterisation of reaction networks that are stochastically complex balanced 
for any choice of rate constants. It is natural to wonder whether a stationary distribution of the form 


(4.3) on some irreducible component P for all choices of rate constants implies something specific about 


the T-system. If for specific form we intend deficiency zero and weakly reversible, this is not the case, 
as this is violated in Example However, in Example the system might be described equivalently 
by means of a weakly reversible deficiency zero system for any irreducible component. The question of 
whether this is always true remains open. We provide here two more examples. 

Example 6. Consider the stochastic mass-action system 


2A ^ 2B 


A + iB ^^Ay B. 


The underlying reaction network is considered in Figure]^ On the irreducible component P = {(1, 5), (3, 3), (5,1)}, 

the Markov chain associated with the system has the same distribution as the Markov chain associated 

with 

2A ^ 2B, 


3K2 


since the transition rates coincide. It is interesting to note that the dynamics of the two systems are 


different when they are deterministically modelled [5]. Due to Theorem 3.3 the latter system is complex 
balanced for any choice of rate constants. Therefore, by Theorem |5.H the stationary distribution on P 


has the form (4.3) on both systems for any choice of rate constants. The same argument does not hold, 


in this case, for the other irreducible components. 
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Example 7. The same phenomenon as in Example is observed in the stochastic mass-action system 

2A^'iA + B A + 3B ^ 2B. 

On the irreducible component T = {{x\,X 2 ) € : Xi> 2, Xi = X 2 }, the Markov chain associated with 

the system has the same distribution as the Markov chain associated with 

2A^2,A + B, 

K 2 

since the transition rates coincide, and the latter network is weakly reversible and has deficiency zero. 


A Preliminary results 


Here we state some preliminary results that will be needed in Appendix |B| 

Lemma A.l. Let Q be a reaction network. // j/i — t 1/2 • • • — >■ is a directed path in the reaction 

graph of Q, and x > yi, then x + yq — yi is accessible from x. 


Proof. First, note that 


q-l 

X + '^{yi+1 -y^) = x + yq- yi. 


2=1 

It is sufficient to note that if x > yi, then for any 1 < j < g — 1, we have 


X + '^{yi+i - yi) = x + yj -yi> yj. 

i=l 


This concludes the proof. □ 

Lemma A.2. Let T be an irreducible component such that Qy has deficiency zero. Then, Qr is weakly 
reversible. In particular, if Q has deficiency zero, Qr has deficiency zero and is weakly reversible for 
every irreducible component F. 

Proof. If T^r is empty then is weakly reversible and there is nothing to prove. Otherwise, if TZr is 
non-empty, let yi —t yi € Hr- By hypothesis, there exists a state x in F with x > yi. This means that 
X + is accessible from x. Moreover, since x belongs to an irreducible component F, we have that 

X is accessible from x + S'® well, which implies that 


9 

x = x + '^i 
i=i 


Vj^vp 


for a certain choice of 


Vi^Vi 


In particular, ~ 0- By the hypothesis of deficiency zero, it 

cause g?, defined in 

— x yj ' Hi ' ’ 

S associated with Qy- Therefore, 


follows that — 0; because (/?, defined in ( |3.3[ ), is an isomorphism between the spaces D and 


9 

~ ^Vj) ~ ^v^y — Oj 
j=i yeCr 

for some integers ay. Since the vectors Cy are linearly independent, Uy = 0 for all y S Cp. Hence, each Cy 
that appears in the sum, must appear at least twice, once with coefficient 1, once with —1. Consequently, 
by iteratively reordering the terms dy.^y'_, the reactions {yj —>■ y'jYj^i form a union of directed closed 
paths in the reaction graph of Q. In particular, the reaction yi —?> y[ is contained in a closed directed 
path of the reaction graph of Qy, and since this is true for every reaction in TZy, Gr is weakly reversible. 
We conclude the proof by Lemma since if Q has deficiency zero, so does every subnetwork of □ 
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Lemma A.3. Let Q be a weakly reversible reaction network, and let T be an irreducible component. 
Then, for any complex y' S Cr we have 

{y &C-. y ^ y' ^11} = {y &Cr-y ^ y & 7^^}, 

{y &C \ y' ^ y & = {y ^Ct- y' ^ y & 7^^}■ 


Proof. One inclusion is trivial, since T^r O TZ. For the other inclusion, fix y' G Cr- Suppose that there 
exists a: S r with x > y'. It follows that any reaction y' —>■ y G TZ is active on F, and therefore is 
contained in T^r- Moreover, since Q is weakly reversible, for any reaction in TZ of the form y —>■ y', there 
exists a directed path in the reaction graph of Q from y' to y. Hence, by Lemma A.l[ x + y — y' is 
accessible from x, which implies that a; + y — y' is in F and that y — > y' is in T^r, since x + y — y' > y. 
Therefore, to conclude the proof it suffices to prove that there exists a; € F with x>y'. 

If it were no x G F with x >y', then no reaction of the form y' ^ y would be in TZ^. Since y' G Cr, 
there exists a reaction of the form y ^ y'. This means that there is i G F, such that x > y. Hence, 
x + y' — y is in F with x + y' — y > y', which concludes the proof. □ 


B Proofs 


B.l Proof of Theorem 13.II 

It is proven in m that if a deterministic reaction system {Q, K) is complex balanced, then Q is weakly 
reversible. By m. we also know that if K is mass-action kinetics, then all positive equilibria are 
complex balanced, and there exists exactly one positive equilibrium in each stoichiometric compatibility 
class, which is locally asymptotically stable. Therefore, to conclude the proof we only need to prove that 
in a complex balanced mass-action system {Q, k), the eventual equilibria on the boundary of K" are also 
complex balanced. 

First of all note that any subsystem {Gl,kl) of {G,k) corresponding to a linkage classes L of G is 
complex balanced. Indeed, the projection of a positive complex balanced equilibrium of {G, k) onto the 
space of the species of L satisfies (3.2) for any complex of Gl, hence it is a positive complex balanced 
equilibrium of {Gl,kl). 

Let c be an equilibrium point on the boundary. Consider a linkage class L of G, and assume that 
cs > 0 for any species S appearing in the linkage class. Then, the projection of c onto the species of L 
is a positive equilibrium of {Gl,kl), and therefore complex balanced. It follows that c satisfies (3.2) for 


any complex of L. Oppositely, assume that there exists a species appearing in the linkage class L, such 
that cs = 0 (this can only happen on a boundary state). Remember that by mass-action kinetics, all 
the rates of reactions whose source complex contains S are zero. In particular, all the rates of reactions 
degrading S are zero. Consider a complex y in L that contains S. By weakly reversibility, there exists 
a reaction y' —>■ y in L. If y' contains S, then Xyi^y{c) = 0. If y' does not contain S, then the reaction 
y' ^ y produces S. Since the rate of all reactions degrading S is zero at c and c is an equilibrium, then 
Xy'^y{c) must be zero as well. By mass-action kinetics, this means that there exists a species S' ^ S 
such that S' appears in y' and cs' = 0. By iteratively applying the same argument with the new species 
S' and by weakly reversibility, we obtain that Xy^yi{c) = 0 for any reaction y —>■ y' in L. It follows that 
c satisfies (3.2) for any complex in L, since the equation reduces to 0 = 0. Equation (3.2) is therefore 
satisfied for any complex of G and c is a complex balanced equilibrium. This concludes the proof. □ 


B.2 Proof of Theorem 13.41 

By [TOl Theorem 6.1.2], if x G Kg is an equilibrium point and y ^ y' G TZ, then suppy C suppx 
only if y —y' is terminal. Moreover, if suppy C suppx, then suppy C suppx for every complex y of 
~ (4);,Cy, 7?.y). 

Now, suppose that K is mass-action kinetics with rate constants k, and that suppy C suppx with 
y y' GTZ (and therefore y ^ y' G TZ*). Consider 

= {y ^ y' & TZ: suppy C suppx}. 

By the first part of the statement, the reaction graph of the subnetwork G = {X,C,TZ) is a union of 
terminal strongly connected components of G, and therefore G is weakly reversible. Moreover, by Lemma 
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3.2 


the deficiency of Q is 0. It is not hard to see that the canonical projection of x onto the spa ce of 


the species A” is a positive equilibrium point of {Q, k), and therefore complex balanced by Theorem 
The proof is concluded by (|3.2[) and by noting that, for any complex y G C- 


3.3 


■-y, 


{y' G C: y ^ y' G TZ} — {y' G Cy : y ^ y' G TZy}, 
{y' G C : y' ^ y G TZ} = {y' G Cy : y' ^ y G TZy}. 


□ 


B.3 Proof of Proposition |4.2| 


If TZr is empty there is nothing to prove. Suppose that this is not the case. Since has deficiency zero, 

y' 


by Lemma A.2 it is weakly reversible. For any y ^ y' G TZr, by definition there exists x S T such that 
X > y-, which in turn implies x y' — y > y' . Therefore, for any directed path in the reaction graph of Q 
that starts with y ^ y' G TZr, all the reactions in the path belong to TZr, by definition of TZr- Since 
is weakly reversible, this can only happen if T^r ^ TZ*, and this proves the first part of the statement. 


To conclude the proof, note that if the deficiency of Q is zero, then by Lemma 3.2 the deficiency of Qr 
is zero as well. □ 


B.4 Proof of Theorem 14.31 

For the first part of the statement, consider a continuous-time Markov chain Cr{t) with state space F x C 
and transition rate from (x, y) to {x + y' — y, y') given by Xy^yr{x) H y ^ y' G TZr, and zero otherwise. 
The master equation for C'r(t) is 


yeCr 


Tt{x - y' + y,y)Xy 


,{x-y' + y)= ^ n{x, y')Xy'. 
yeCr 


,{x) My'GC,xGT, 


with the convention that Xy^y'ix) = 0 if y —> y' ^ T^r- By Definition a stationary distribution for 
Cr{t) exists and it is of the form TT{x,y) = Mtt(x), for a suitable normalising constant M. Since 7r(x) 
is positive for any a; S F (because it is a stationary distribution on an irreducible component), then by 
standard Markov chain theory, we have that for any two states (a;i,yi), {x 2 ,y 2 ) S F x C, if (a; 2 , j/ 2 ) is 
accessible from (a;i,yi), then (xi,yi) is accessible from ( 2 : 2 , ^ 2 )- Fix y ^ y' G TZr and a; G F with x >y. 
Then, a directed path from {x + y' — y,y') to {x,y) exists in the graph associated with Cr{t). The second 
components of the form y of the states in the path, by construction, determine a directed path in the 
reaction graph of Qr from y' to y. Hence, any reaction y ^ y' G TZr is contained in a closed directed 
path, which means that is weakly reversible. 

Assume now that K is mass-action kinetics with rate constants k and that c is a positive complex 


balanced equilibrium of {Q, k). Then, by Theorem 4.1 there exists a (unique) stationary distribution on 


F of the form (4.3). If a species Sj is not in Xr, then the value of Xj is constant for any x G F, and (4.5) 


can be obtained from (4.3) by modifying the normalising constant. 


By Theorem |3.1| and Lemma [A.3[ we have that 


E- 

yeCr 


iy^y' 


= E' 

yeCr 


'^y'—^y '^y ^ Fp 5 


with 


ly^V' 


= 0 if y —y' ^ T^r- Therefore, for any y' G Cp and x G F, 


-i_y 


x+y-y' 


(x - y')\ 




!/6Cr 


^ ^ '' y&Cr 


which leads to (4.4), since tt is of the form (4.3). 

To prove the converse we first introduce a new stochastic mass-action system (^r, kp), which is given 
by the reactions of the form 


y + Sy ^ y' + Sy! with y y' G TZr, 
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where Sy are fictitious species in one to one correspondence with the complexes Cr- The rate constant 
of the reaction y + Sy ^ y' + Sy' is given by Ky^y' ■ It is not difficult to see that the sum of the fictitious 
species is conserved for any possible trajectory. Moreover, since any directed path j/i —?> ?/2 2/9 in 

the reaction graph of Q corresponds to a directed path yi + Sy^ —>■ 2/2 + Sy^ —)■ .. .yq + Sy^ in the reaction 
graph of Qy y we have that Qy is weakly reversible by the first part of the proof. 

Consider the set 

T = {{x,x) G N" X N"" : a; G r,||i||i = 1}. 

Every state in T is of the form {x,Sy) G where a; G T and Sy is considered as the vector in N™ 

with entry 1 in the position corresponding to the species Sy and 0 otherwise. Since T is an irreducible 
component of Q and the sum of the fictitious species is conserved, no state outside T is accessible from 
any state in T, according to Qy- Moreover, the master equation on T can be written as 


/ '1 c \ ( 2 ^ ~ 2/^ + 2/)! T| 

7r(a; y +y,Sy)Ky^y> _ y/y^ 

x\ 

= 'y^Tr{x,Sy')Ky’^y, V 2 / GC,XGT. (B.l) 

yeCr ^ ^ '' 


V&Cr 


If we choose %(x,x) = Mt:(x) for some positive constant M, then the master equation (B.l) is satisfied 
due to Definition]^ Therefore, if M is chosen as a suitable normalising constant, tt{x,z) = Mtt{x) is a 
stationary distribution on T. 

Consider the linear homomorphism if as defined in (3.3), for the reaction network Qy- Let | • | denote 
the cardinality of a set, and note that |Cr| = |Cr| = wr- For any vector Cy of the basis of K'"'", we have 
ip{ey) = {y,Sy). Since the vectors (y^Sy) with y G Cy are linear independent, ip is an isomorphism and 
the deficiency of is 0 . 


Since is a deficiency zero weakly reversible reaction network, it foll ows from Theorem 


4.1 


3.3 


that 


we have that tt has the 


the mass-action system (^r, k) is complex balanced. Therefore, by Theorem 
form 

Tt{x,x) = 

^ x\ x\ 

for a positive complex balanced equilibrium (c, c), on any irreducible component T contained in T. Since 
7r(x,x) = M7r(x) does not depend on x, we have 


7r(x,x) =M^ — , 


for any (x,x) G T. 

Fix a complex y' G Cy- Since is weakly reversible, there exists a reaction y' y that is active on 
F. Fix a: G F such that x > y' . Then for any y ^ y' G TZy we have x — y' + y>y- If we plug the formula 


for 7r{x,x) in (B.l) for our choice of x and y' , we obtain 


f-x-y'+y 

yeCr 


{x-y' + y)\ 




{x-y' + y)\ 
{x-y')'. 


yeCr 


xl 


' {x - y’)'.' 


which leads to 

yeCr yGCr 

The proof is concluded by the fact that the above holds for any fixed y' G Cy, which means that c is a 
positive complex balanced equilibrium of (^r, ^r). □ 


B.5 Proof of Theorem 14.61 


By Lemma A.2 Qy is weakly reversible. Moreover, for y ^ y' G TI-y, 'ii x > y then x + y' — y > y'- 
This implies that for any directed path in the reaction graph of Q that starts with y ^ y' G T^r, all the 
reactions in the path belong to T^r, by definition of T^r- Since is weakly reversible, every directed 
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path in the reaction graph of G that starts with y ^ y' G TZr is contained in a closed directed path. This 
implies that TZr C TZ*, and proves the first part of the statement. 

Now assume that K is mass-action kinetics with rate constants k. If the deficiency of G is zero, then 
by Lemma |3.2| the deficiency of the terminal network is zero as well. Moreover, G* is weakly reversible 
by definition, thus by Theorem 3.3 {G*, k*) is complex balanced for any choice of rate constants k* . 

Let X{t) be the stochastic process associated with {G, k). By the first part of the statement, on T only 
terminal reactions take place and these involve a subset of the species only. Without loss of generality, we 
can assume that X* is constituted by the first n* species of X. Therefore, T is of the form T* x {u}, with 
r* C 7^" and v € M"”" . Moreover, we have that on T*, the projection X*(t) = {Xi(t),..., Xn*(t)) 
is distributed as the process associated with {G*,k*), for which T* is an irreducible component. Let c 
be a positive complex balanced equilibrium for {G*,^*). Hence, by Theorem 4.1 or Corollary |4.4| 


stationary distribution of the process X(t) = {X*{t),v) on T is of the form (4.6). 


the 

□ 


B.6 Proof of Theorem 14.71 

For the first part, we prove that if an irreducible component T is positive, then G is weakly reversible. 
This simply follows from Lemma |A.2| indeed, by the lemma, is weakly reversible and since F is 
positive, Gr = G- 

To prove the second part, we have to show that a weakly reversible reaction network is essential, 
and this is done in [22) . Moreover, a deficiency zero weakly reversible mass-action system is complex 
balanced, and the proof is concluded by Theorem |4.1| or Corollary |4.4| □ 


C Calculations for Examples and 

In Example we claim that the stationary distribution on the irreducible component Fg = {a; S 
: rci -I- 2:2 = 0} has the form 


TTe{xi,X2) = Me — j —- for (si, 2 : 2 ) G Fg. 
a;i!a;2! 

To prove this, it is sufficient to show that Trg satisfies the master equation for every point (xi, 2 : 2 ) of Fg. 
The master equation on (xi,X 2 ) is given by 


pTTg{xi -I- I,X 2 - 1)(6» - I)(xi -I- I)1{2,2>1} -I- pTT 0 {xi - 2, X 2 + 2) 1 > 2 } 

X 2 ■ 

/ (X 2 )’ 

= p7re(xi,X2) ((61 - -k _ 2)1 

By plugging in the formula for irg and after dividing by p and Mg we obtain 

— ]—Ax2{9 - 1) -k xi(xi - I)] = —^[a;i(6' - 1) -k X 2 {x 2 - I)]. 

Xi!x2! xi!x2! 

If we multiply by xi!x 2 ! and substitute 0 = xi -k X 2 , it follows that 

X 2 {xi -k X2 — I) -k Xi(xi — 1) = Xi(xi -k X2 — I) -k X 2 {X 2 — 1), 


that is 

X1X2 + x\ — X2 + x\ — Xi = x\ + X1X2 — Xi + x\ — X2, 

which always holds true because the terms cancel each other. 

In Example]^ we change the notation to Fg = {x G : xi -k X 2 = 0 -k 1}. Then we claim that the 

stationary distributions on the irreducible components Fg^ and Tg^ are irg^ and Trg^, respectively, where 

as before ^ 

■ne{xi,X 2 ) = Mg —,— , for (xi,X 2 ) G Fg. 
xi!x2 ! 


18 




















We prove that irg^ is the stationary distribution on rg^. The case with 62 is analogue. We prove the 
result by consider the master equation for Trg^ on a point (xi,X 2 ) G Tgj^, which is as following: 


pireA^i + l,X2- l)0i02ixi + l)l{2;2>i} 


+ — 2 ,X 2 + 2 ){ 0 i + 02 — 1 )^^-f^^{a:i> 2 } 

(xi + 2 )! 

+ pT^SiiXl + 2 , X 2 — —— 1)! 

I _ / ,0 o^ (^1 + - 2 )! ^ 

+ P7rg^(a:i + 2 ,X 2 — 2 ) — xi\{x 2 — 3 )!-^{a;i>o,a;2>3} 

0 : 2 ! 


= p'n9^{xi,X2)0i02Xil{xi>i} + p7rg^(a;i,X2)(6»i +02- 1) ^^ _'2)! 

X\ ! X\ \x 2 • 


As we did for the previous calculations, we plug in the expression for TTg^, then divide by Mg^, p and 
multiply by Xi\x 2 \- We obtain 

0102 X2 + (01 + 02 - l)a;i(a;i - 1) + XiX2{X2 - 1) + X2ix2 - I){x2 - 2) 

= 01^22:1 + {01 + 02 — ^)x2{x2 ~ 1) + 2:i(2:i — l)(2:i ~ 2) + Xi{xi — l)a:2. 

Finally, by substituting 0i with cci + 0:2 — 1 and by performing the calculations, we obtain 0 = 0, which 
means that the above equation is satisfied. 


References 

[1] David F. Anderson, Gheorghe Craciun, Manoj Gopalkrishnan, and Garsten Wide, 
Lyapunov functions, stationary distributions, and non-equilibrium potential for chemical reaction 
networks. arXiv: 1410.4820, 2014. 

[2] David F. Anderson, Gheorghe Graciun, and Thomas G. Kurtz, Product-form stationary 
distributions for deficiency zero chemical reaction networks, Bull. Math. Biol., 72 (2010), pp. 1947- 
1970. 

[3] David F. Anderson, German A. Enciso, and Matthew D. Johnston, Stochastic analysis 
of biochemical reaction networks with absolute concentration robustness, J. of the Royal Society 
Interface, 11 (2014), p. 20130943. 

[4] David F. Anderson and Thomas G. Kurtz, Continuous time Markov chain models for chem¬ 
ical reaction networks, in Design and Analysis of Biomolecular Circuits: Engineering Approaches 
to Systems and Synthetic Biology, Heinz Koeppl, Douglas Densmore, Gianluca Setti, and Mario 
Di Bernardo, eds., Springer, 2011, pp. 3-42. 

[5] David F. Anderson and Thomas G. Kurtz, Stochastic analysis of biochemical systems. 
Springer, 2015. 

[6] Karen Ball, Thomas G. Kurtz, Lea Popovic, and Greg Rempala, Asymptotic analysis of 
multiscale approximations to reaction networks, Ann. Appl. Probab., 16 (2006), pp. 1925-1961. 

[7] L. Boltzmann, Neuer beweis zweier sdtze iiber das wdrmegleichgewicht unter mehratomigen gas- 
molekiilen, Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften in Wien, 95 (1887), 
pp. 153-164. 

[8] Gheorghe Craciun and Casian Pantea, Identifiability of chemical reaction networks, J. Math. 
Chem., 44 (2008), pp. 244-259. 


19 








[9] Peter Erdi and Janos Toth, Mathematical models of chemical reactions: theory and applications 
of deterministic and stochastic models, Manchester University Press, 1989. 

[10] Martin Feinberg, Chemical reaction network structure and the stability of complex isothermal 
reactors — I. The deficiency zero and deficiency one theorems, Chemical Engineering Science, 42 
(1987), pp. 2229-2268. 

[11] Ankit Gupta and Mustafa Khammash, Determining the long-term behavior of cell populations: 
A new procedure for detecting ergodicity in large stochastic reaction networks, in Proceedings of the 
19th IFAC World Congress, 2014. 

[12] K. L. Hey, H. Momiji, K. Featherstone, J. R. E. Davis, M. R. H. White, and D. A. 
Rand, A stochastic transcriptional switch model for single cell imaging data. Biostatistics, in press 
(2015). 

[13] Fritz Horn and Roy Jackson, General mass action kinetics. Arch. Ration. Mech. Anal., 47 
(1972), pp. 81-116. 

[14] Piers J. Ingram, Michael P.H. Stumpf, and Jaroslav Stark, Nonidentifiability of the source 
of intrinsic noise in gene expression from single-burst data, PLoS Computational Biology, 4 (2008), 
p. el000192. 

[15] J. R. Jackson, Networks of waiting lines, Oper. Res., 5 (1957), pp. 518-521. 

[16] Hye-Won Kang and Thomas G. Kurtz, Separation of time-scales and model reduction for 
stochastic reaction networks, Ann. Appl. Probab., 23 (2013), pp. 529-583. 

[17] Frank P. Kelly, Reversibility and stochastic networks, John Wiley & Sons, Inc., 1979. 

[18] Thomas G. Kurtz, The relationship between stochastic and deterministic models for chemical 
reactions, J. Chem. Phys., 57 (1972), pp. 2976-2978. 

[19] -, Strong approximation theorems for density dependent markov chains. Stochastic Process. 

Appl., 6 (1978), pp. 223-240. 

[20] Jean Mairesse and Hoang-Thach Nguyen, Deficiency zero petri nets and product form, in 
Applications and Theory of Petri Nets, Springer, 2009, pp. 103-122. 

[21] Andrea Marin, Simonetta Balsamo, and Peter G. Harrison, Analysis of stochastic Petri 
nets with signals. Performance Evaluation, 69 (2012), pp. 551-572. 

[22] Loi'c Pauleve, Gheorghe Graciun, and Heinz Koeppl, Dynamical properties of discrete re¬ 
action networks, J. Math. Biol., 69 (2014), pp. 55-72. 

[23] Peter Whittle, Systems in stochastic equilibrium, John Wiley & Sons, Inc., 1986. 

[24] C. Zechner, M. Unger, S. Pelet, M. Peter, and H. Koeppl, Scalable inference of heteroge¬ 
neous reaction kinetics from pooled single-cell recordings, Nat. Methods, 11 (2014), pp. 197-202. 


20 



